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P3 ' A numerical method is devised for study of stochastic partial differential equations 

CJ ■ (SPDEs) describing directed percolation, the contact process, and other models with 

a continuous transition to an absorbing state. Owing to the heightened sensitivity to 
;_j ■ fluctuations attending multiplicative noise in the vicinity of an absorbing state, a useful 



C^ 



method requires discretization of the field variable as well as of space and time. When 
applied to the field theory for directed percolation in 1+1 dimensions, the method 
yields critical exponents which compare well against accepted values. 



PACS numbers: 05.50. +q, 02.50.-r, 05.70.Ln 



1 Introduction 

The study of critical phenomena in simple nonequilibrium lattice models has reached the 
stage where many transitions can be assigned to one of a small set of universality classes. 
For continuous transitions into an absorbing state, a very high degree of universality has 
been found, with many examples supporting the prediction [|l], ^, ^ that such transitions 
belong generically to the class of directed percolation (DP). Examples include the basic 
contact process and its variants @, |^, H, 0, H], surface reaction models |P, |^, 0], branching 



and annihilating random walks with odd parity [11, 12, 13], assorted multiparticle processes 



P^ P^ , [T3| , P^ , and even models with multiple absorbing configurations |TB|, |T^, ^ ^, |2^ . 
Each is an interacting particle system characterized by rules for elementary processes such as 
creation, annihilation, and diffusion. Looking at the rules, there is little to tell us what sort 
of critical behavior to expect, nor why it is universal. Understanding of universality emerges 
instead from the study of coarse-grained formulations which capture the large-scale features 
essential to critical behavior. In such field theories the microscopic picture of particles on a 
lattice is replaced by a set of densities which evolve via stochastic partial differential equations 
(SPDEs). At this level, renormalization group methods may be applied, fl], |2^, |2^, |2^, |26[| . 



A basis for universality appears if one can show that the continuum descriptions for various 
models differ only by irrelevant terms. At present, however, there are many more models 
known (on the basis of numerical work - simulation and/or series analysis), to have DP 
critical behavior, than have been studied using field theory. Useful continuum descriptions 
of multiparticle processes, for example, have yet to be devised 0, |??[] . 

It is of interest, therefore, to study SPDEs for nonequilibrium systems, and to compare 
their behavior with the lattice models they are supposed to describe. But solving a non- 
linear SPDE is not generally feasible by analytic means, and so numerical methods must 
be sought 1^. Numerical integration has been applied to several SPDEs, for example the 
time-dependent Ginzburg-Landau equation describing phase separation BT], B^, P3[, and the 



Kardar-Parisi- Zhang equation pj, 35, p&. In problems with an absorbing state, however. 



the usual approach does not yield useful results. A method for dealing with such systems is 
proposed in the present work, and is used to study the field theory for the contact process. 

The outline is as follows. In Section 2 I describe the original model and the corresponding 
SPDE. The integration scheme is introduced in Section 3, and results are presented in Section 
4. A discussion and summary follow in Section 5. 



2 Lattice Model and Field Theory 

In the contact process (CP) ||^, each site of the d-dimensional cubic lattice, Z'^, is either 
vacant, or is occupied by a particle. The transition rules are easily stated: a vacant site with 
n occupied nearest neighbors becomes occupied at rate Xn/2d, and particles disappear at 
unit rate, independent of their surroundings. Evidently the vacuum is absorbing; the active 
phase, characterized by a nonzero stationary particle density, p, exists only for sufficiently 
large A (and only, strictly speaking, in the infinite- volume limit). There is a continuous 
transition from the vacuum to the active phase at a critical value Ac [0 . (In one dimension 



Ac — 3.2978 [0.) The transition belongs to the universality class of directed percolation. 



(Note that the d-dimensional CP corresponds to directed percolation in d+1 dimensions.) 

Janssen [0] proposed a continuum description of the CP and allied models: 
ap(x, t) 



dt 



ap{x, t)-hp^ -cp-^ + --- + DV'p + r/(x, t). (1) 



p(x, t) > is the coarse-grained particle density; the ellipsis represents terms of higher order 
in p. ri{'x,t) is a Gaussian noise, which respects the absorbing state (p = 0) by virtue of the 
covariance: 



r7(x, t)r]{x', t') oc p(x, t)(5(x - ^)5{t - t'). (2) 



This form can be justified by coarse graining the CP, in the hmit of large bin size. Let 
Tii be the number of particles in bin i, and An\ the change in this number during a brief 
interval. The latter has expectation Arii oc arii + 0{nf)^ (with a oc A — 1), and under the 



customary assumption of Poissonian statistics for reaction systems, its variance equals Ani. 
For sufficiently large bins we may approximate Ariiby a Gaussian. Thus, since reactions in 
different bins are uncorrelated, coarse-graining the original model leads to a stochastic field 
theory with Gaussian noise whose autocorrelation is proportional to the local density. (There 
is also noise due to the fiuctuating diffusive current. But diffusive noise does not affect the 
critical behavior in the present case, and so I shall ignore it, in the interest of simplicity.) 
Since Eq (|l|) involves multiplicative noise, one must decide upon an interpretation ||29|. As 
shown in the following section, the Ito interpretation of Eq (|I|) is demanded by physical 
considerations. 

In mean-field approximation (the spatially uniform, noise- free version of Eq (|l|)), the 
vacuum becomes unstable when a = 0, and for a, 6 > there is an active state. When 
fluctuations are taken into account, the critical point shifts to Oc > 0, and the critical 
behavior is nonclassical. For example, the stationary density in the CP scales as poc {a — acY , 
with (3 ~ 0.277 in one dimension. (In mean-field theory, (3 = 1.) Field-theoretic analysis 
ll], Q reveals that the cubic and higher-order terms are irrelevant to critical behavior, so 
long as 6 > 0. (Such terms are therefore ignored in what follows.) The situation is analogous 
to that in equilibrium critical phenomena, where the Ising universality class is generic for 
models with a scalar order parameter and short-range interactions. 

Without noise, Eq (|T]) is a reaction-diffusion equation, which exhibits a mean-field critical 
point. It is perhaps surprising that driving a react ion- diffusion equation with multiplicative 
noise leads to the proper exponents. Of course the condition expressed in Eq (^ is crucial 
in this regard. On the other hand, it is not clear whether adding a properly scaled noise to 
the reaction diffusion equation always yields a useful field theory [H, ^]. 

Further unanswered questions are whether solutions to Eq (|T]) exist, and if so, whether 



they reproduce the phenomenology of the lattice models they are supposed to describe. 
(For example, Can the field actually fiuctuate into the vacuum?) Such issues never arise in 
renormalization group analyses, where the SPDE merely serves as a basis for perturbation 
theory, which proceeds by expanding the formal solution. Since the exponents emerging from 
the e-expansion analysis of Eq (|I]) are in good agreement with series and simulation results, 
there is no reason to doubt its validity in this context. The present work is concerned with 
nonperturbative (numerical) solutions to a discretized version of the SPDE. 



3 Numerical Method 

Can Eq ([1|) be integrated numerically? To begin, we discretize space, obtaining a set of 
Langevin equations which in one dimension take the form 

^^ = ap{^, t) - bp' + DV'p + r/(z, t), (3) 

where z is a site index, (for convenience we assume a spacing Ax = 1 in the discretization), 
and V^p{i, t) = p{i + 1, t) + p{i — 1, t) — 2p{i, t) is the lattice Laplacian operator. The noise 



term satisfies 77(2, t)ri{j, t') = Tp{i, t)6ij6{t — t'). iV is related to the growth rate A in the CP. 
We may regard it as constant over the range of parameter values of interest here, and set 
r = 1 from here on.) Applying the Cauchy-Euler method to these equations [^, we find 



p(z, t + At) - p(z, t) = [ap{i, t) - hp{i, ty + DVy{i, t)]At + ^ p{i, t)AtY{i, t). (4) 

where the Y{i,t) are independent Gaussians with zero mean and unit variance. Eq (^ is 
similar to the set of stochastic differential equations (SDEs) derived in Ref.[^] from dis- 
cretization of a SPDE. The latter scheme, however, is not useful here, due to the dominance 
of the noise in the vicinity of the critical point. Indeed, once we discretize time there is noth- 



ing to prevent p{i, t) becoming negative. We might attempt to remedy this by stipulating 
that whenever integration yields p{i,t) < 0, the density at that site be set to zero. But this 
artifice is not without drawbacks. We are interested in the critical region, where p is small. If 
the typical magnitude of the first term on the r.h.s. is e, the noise term is of order y/e, and so 
overwhelms the deterministic part of the evolution. In the original equation, the cumulative 
effect of the systematic term is not obliterated by the noise, which has zero mean. But this 
is no longer so if we truncate the noise in an unsymmetric manner. In simulations using the 
max[p, 0] rule, the process never reaches the vacuum (even for a < 0). Regions of density 
zero are rapidly repopulated by nearby active sites. It appears, then, that straightforward 
numerical integration of Eq (^ is not useful. 

Consider, for the moment, the Ito SDE (a Malthus-Verhulst process), with the same local 
terms as Eq ([l|): 

dp = [ap — hp^]dt + y/pdC,it), (5) 

with dC,{ty = dt. Eq (^ describes Brownian motion in a potential which grows oc p^ for large 
p, and has a minimum (for a, 6 > 0), at p = a/h. There is an absorbing boundary at p = 0, 
corresponding to the vacuum in the CP. Now suppose we had interpreted Eqs (|l|) and (^) 
as Stratonovich equations. The Ito SDE corresponding to the Stratonovich interpretation of 



Eq d) is [H 



dp = [ap - hp^ + -Adt + ^/pd^{t). (6) 

It includes a constant source term, so that p = is no longer absorbing. Clearly this is not 
the problem we began with! Hence Eqs (|I]) and (|^) should be taken in the Ito sense. 

Integrating Eq(H), we have: 

Ap = [ap - bp^]At + ^AW, (7) 



where AW = yAtY, and Y is Gaussian with zero mean and unit variance. Now to prevent 
p + Ap going negative, I propose to discretize the density by setting p = npmin, in > 0), and 
at the same time to truncate Y symmetrically by restricting its magnitude so: \Y\ < Y^ax- 



We require l^aa;V At < ^ pmm to avoid negative densities. This can be achieved in a variety 
of ways, for example by setting 



In At 

Y = - 

^ max 

and 



Y — 

-^ max Q ; 



Eqs (D and (^ represent but one of an infinity of choices. Determining which is optimal 
for a specific problem is left as a subject for future work. I take pmin oc At in hopes of 
minimizing the effect of a discretized density. The relatively slow growth of Y^ax poses no 
essential difficulty. (Note that for At = 10~^ we have Y^ax = 3.07 — about three standard 
deviations.) Indeed, all noise distributions having zero mean and finite variance should 
yield qualitatively similar behavior. If one were interested solely in universal properties, the 
Gaussian could be replaced with a uniform distribution, in the interest of computational 
efficiency. 

Having discretized p, we can define an integer process by exploiting the invariance of Eq 
(0) under the rescaling: 



(10) 
(11) 
(12) 

If we choose a = Pmin, then p' is restricted to integers > 0. Discretization (in time) of 



b- 


^b' = ah 


p- 


, P 


e- 


a 



Eq (1^) leads to a noise term Yy/pAt; in the rescaled equation it becomes YJp'At/a. (Y is 
a zero-mean, unit- variance Gaussian, truncated as described above.) 

We now have a discretized version of Eq (|^), in which positivity and zero-mean noise are 
ensured at the cost of a "quantized" density. Since p' can change only by integer steps, it is 



likely (especially for small p'), that many increments of the density will be rejected, for being 
of less than unit magnitude. It therefore seems advisable to introduce a continuous variable 
ip which accumulates the increments in density at each time step. Whenever {ipl > 1, the 
integer part is transferred to p'. 

In summary, the numerical scheme for Eq (^) is as follows. At each time step tp — > ip+Atp, 
where 

All; = {ap' - b'p'^)At + Y^p'At/a, (13) 

and 

P'-P' + M, (14) 

^^^-[ip], (15) 

where square brackets denote the integer part. (Initially, tp is zero.) Eqs (|TB|) - (|T^ may 
be viewed as a Malthus-Verhulst process in which the population change. An = Ap', is 
approximated by a suitably truncated Gaussian random variable. Simulations show the 
population fluctuating around a quasi-steady value, Pgs ~ a/b, and eventually becoming 
trapped at zero (see the inset of Fig 1). In Fig 1 the mean density for a sample of 10^ trials 
is plotted for several time increments. (The model parameters are a = 1.5, b = 1; p = 1.6 
initially.) The density decays exponentially, with relaxation times 12.9, 9.7, and 9.8 for time 
increments 10^^, 2 x 10~^, and 10~^, respectively. This is in good agreement with the mean 
first passage time, 10.3, for hitting p = 0. (The latter is obtained from the Fokker-Planck 
equation corresponding to Eq (^) pp|.) 



Our treatment of the SPDE, Eq dTj), closely parallels that of Eq (^). Discretization and 
rescaling of Eq(H) yields a set of diffusively coupled Malthus-Verhulst processes: 



Aij, = {ap[ - b'pf + DV'p[)At + Y,^p'At/a, (16) 

and 

p--p- + m, (17) 

^i^V'^-M- (18) 
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(The rescaling of Eqs (|iy) - (p!2| ) does not affect the diffusion coefficient.) We have con- 
verted the original SPDE into a lattice of discrete stochastic processes which approaches the 
continuum model as At and Ax -^ 0. 

A final technical point is that when integrating the coupled equations, the three parts 
of the evolution — deterministic on-site contributions, the noise term, and diffusion — are 
implemented separately, in turn (over the entire lattice), at each step. (Transfer from ip to p' 
is made following each of the three sub-steps.) Unphysical events such as a site becomming 
empty, and subsequently acting as a source for its neighbor, are eliminated by this measure. 



4 Results 

I applied the scheme detailed above to systems of several hundred to several thousand sites, 
in order to study the critical behavior of Eq (|I|). To begin, I describe results for the steady- 
state, obtained in a series of long runs (of duration tf ^ 10^), on lattices of 500 — 2000 
sites, using a time step At = 10~^. The coefficient b was set to unity; diffusion rates D = 1 
and 10 were considered. The stationary density p(a) (expressed in its original units, prior 
to rescaling), is shown in Fig 2. The data suggest a continuous transition to the vacuum at 
a critical value Qc (— 0.77, 0.36, for D = 1, and 10, respectively). To estimate the order- 
parameter exponent (3, one must estimate the critical value Oc, and then plot the density 
versus A = a — Oc on log scales. For Z) = 1, a reasonably linear plot (for small A) is obtained 
when we choose Oc = 0.769 (see Fig 3). A least-squares linear fit to the five points nearest 
Cc yields a slope f3 = 0.295, in fair agreement with (3 = 0.277 for the one-dimensional CP. 
As is often the case, the slope depends quite sensitively upon one's estimate of the critical 
point, and so this analysis is not very precise. (Taking Oc = 0.765, for example, one finds 
/3 ~ 0.39.) The slope (~ 0.45) obtained from the D = 10 data suggests a more mean-field 
like behavior for faster diffusion. Indeed, ac appears to shift towards its mean-field value, 0, 
with increasing D. For larger D a crossover from mean-field to DP-like behavior presumably 



occurs very near the critical point. 

In order to derive quantitative results on critical behavior, I turn to the "time- dependent" 
method 0, ^ , in which one studies the dynamics of spreading from a distribution localized 



near x = 0. The quantities of interest are the survival probability, P{t), mean total density 
n(t), and mean-square spread R'^it), for a large sample of independent trials, all with the 
same initial condition. P{t) denotes the probability of not being in the vacuum state at 
time t, n{t) is the sum of the site densities (averaged over all trials, including those which 
have reached the vacuum by time t), and -R^(t) = J2j j'^p{j,'t)/n{t). In a subcritical system, 
(a < Oc), we expect P and n to decay exponentially, whilst -R^(t) oc t. For a > ac, P 
approaches a nonzero limiting value, n{t) ~ f^ (in d dimensions), and -R^(t) ~ t^, as a 
fraction of trials survive indefinitely and spread at finite speed into the surrounding vacuum. 
At the critical point there is no characteristic time-scale for relaxation, and the evolution is 
characterized by nontrivial power laws: 

P{t) oc t~\ (19) 

n{t) oc f?, (20) 

and 

R\t) oc t\ (21) 

I studied spreading in simulations beginning with a localized density (typically, p{i, 0) = 
Pmin over the 10 - 20 sites nearest the origin, and p{i,0) = elsewhere), ioi D = b = 1. 
Time increments of 10~^ and 10^^ were employed, as At = 10~^ resulted in excessive run 
times. The trials (on the order of 10'^ at each a- value of interest), ran to a maximum time 
of 4000 (1000 in the studies employing At = 10~^), and were performed on lattices of 500 - 
900 sites, sufficiently large that the active region did not reach the boundaries. 

Using the criterion of asymptotic power laws at the critical point, I estimate ac = 
0.7210(5) for At = 10-^ and Oc = 0.568(1) for At = 10"^. A plausible explanation for the in- 
crease in Oc, as At is reduced, is that as pmin becomes smaller, and the truncation of the noise 

10 



less severe, fluctuations to zero density occur more readily. When plotted versus At~^/^, the 
data for Oc suggest a finite limiting value ac(0) ~ 0.8, with ac(At) ^ ac(0) + const, x y/Kt. 
The present very limited data (three points, for D = b = 1), are of course insufficient to 
permit a firm conclusion in this regard. 

Figs 4 and 5 show the evolution of the survival probability, mean population, and root- 



mean-square spread {xrms = \J R'^{i))i ^oi At = lO^"' and 10^^, respectively. In the latter 
case data from slightly off-critical values are also plotted, which show some curvature. From 
least-squares fits to the asymptotic linear region in these log-log plots, I derive the exponent 
estimates listed in Table 1. Uncertainties, given by the figures in parentheses, are subjective 
estimates based on the spread of exponent values found in simulations with a ^ a^. While 
not of high precision, the exponents found here compare rather well against the known DP 
values. Derivation of truly precise results will require longer runs and larger samples, so that 
ttc can be fixed more reliably, and short-time corrections to scaling can be eliminated, by 
means of a local-slope analysis. 

Fig 6 shows the evolution of the mean density profile in the critical system (At = 10"^, 
a = 0.721). Following an early build-up in the central region, the profile broadens and 
becomes more sparse. An interesting aspect of the late-time profile, which merits further 
study, is its large spread, compared to a Gaussian distribution. (That is, x'^ > 3x'^.) 



5 Discussion 

We have seen that some care is required in integrating a field theory with multiplicative 
noise and an absorbing state. To avoid negative densities and complete dominance of noise, 
the SPDE must be regularized in some fashion. The present work shows discretization of 
the field variable to be a suitable method for tempering the equation. Similar conclusions 
apply to the associated SDK. In fact, the method devised here yields an accurate relaxation 
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time for the latter problem. 

Despite discretization of the density, the present approach retains the essential features 
of a continuum description. The density approximates a continuous variable, and spatial 
coupling accurs solely through diffusion. Moreover, creation and annihilation are here ex- 
pressed in a naive mean-field-like manner (they are represented, that is, by terms oc p and 
p^). Nontrivial critical behavior arises by virtue of properly scaled, multiplicative noise. 
The exponent values derived from the discretized SPDE are in rather good agreement with 
accepted values for DP in 1+1 dimensions, arguing for the reliability of the method. The nu- 
merical scheme proposed in this work may therefore be of value in testing candidate theories 
for models which have so far resisted analysis in continuum representation. 
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Table 1. Critical exponents from numerical integration of the SPDE. 



6 


V 


Z/ 


directed percolation 


0.1597(3)'^ 


0.317(2)'' 


1.272(7)'' 


present work: At = 10^^ 


0.15(1) 


0.28(1) 


1.18(2) 


present work: At = 10^^ 


0.159(6) 


0.326(10) 


1.23(2) 



"Ref.f38l "Ref. 41 
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Figure Captions 

Fig. 1. Evolution of tlie mean density in the discretized SDK, Eq (J^), for a = 1.5, 6 = 1. 
Squares: At = 10"^; +: At = 2 x 10^^; o: At = 10""^. The inset shows a typical trial 
(At= 10-3). 

Fig. 2. Steady state density versus a in simulations of the discretized SPDE, Eq (15), with 
6 = 1, and At = 10~^. Solid squares: D = 1; open squares, D = 10. 

Fig. 3. The data of Fig. 2 {D = 1) plotted versus A = a — Oc, assuming Oc = 0.769. The 
straight line, fitted to the five points nearest Oc, has slope 0.295. 

Fig. 4. Time-dependence of the suvival probability, P, total density, n, and mean-square 
spread, Xrms, for a = ac = 0.721, b = D = 1, and At = 10^'^. The straight lines are 
least-squares fits (for slopes see Table 1). 

Fig. 5. Same as Fig. 4, but for At = 10^^. Dots: a = 0.565; open squares: a = 0.568; 
circles: a = 0.570. 

Fig. 6. Average density profiles for D = B = 1 and a = ac = 0.721, At = 10"^. From 
narrowest to most broad: t = 0, 500, 1000, 4000. 
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